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Abstract 

We present a new method for realizing the adiabatic connection approach in density func- 
tional theory, which is based on combining accurate variational quantum Monte Carlo cal- 
culations with a constrained optimization of the ground state many-body wavefunction for 
different values of the Coulomb coupling constant. We use the method to study an elec- 
tron gas in the presence of a cosine-wave potential. For this system we present results for 
the exchange-correlation hole and exchange-correlation energy density, and compare our 
findings with those from the local density approximation and generalized gradient approx- 
imation. 
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1 Introduction 



Density functional theory (DFT) (|l], 0) is the main computational tool for 
the treatment of many-body effects in solid state electronic structure calcula- 
tions and is now widely used to determine ground-state properties of atoms 
and molecules (^. In the Kohn-Sham formulation of DFT (Q) the problem of 
finding the ground state energy and density of an interacting N-electron sys- 
tem is tranformed into an equivalent problem involving non-interacting elec- 
trons. The central quantity in this formulation is the exchange-correlation 
energy E^c, which is a universal functional of the electron density n{T). The 
exchange-correlation energy functional, a complicated many-body object, is 
the big unknown of the theory and the core problem in the density functional 
approach is to find accurate approximations for E^c- The most frequently 
used approximations to date are the local density approximation (LDA) (||) 
and various generalized gradient approximations (GGA) (^, |^, |^). 

An entirely different approach to the ground-state many-body problem is 
quantum Monte Carlo (QMC) (^. QMC calculations are computationally 
more demanding than density functional calculations. However, unlike the 
density functional approach, in which the ground-state density is the basic 
variable, quantum Monte Carlo methods focus on sampling the full ground- 
state many-body wavefunction of the system under consideration and hence 
yield a more detailed description of many-body effects. Quantum Monte Carlo 
calculations can therefore be used to investigate density functional theory from 
"outside" and to test the performance of approximations to E^c- In the last 
few years a number of quantum Monte Carlo investigations of DFT have been 
reported for atoms and molecules (^, |^), model solids (|TU|, [TI|) and silicon (|T^). 
Most of these investigations focused on extracting the exchange-correlation 
potential and components of exchange-correlation energy from accurate elec- 
tron densities obtained from Monte Carlo calculations. Except for a very 
recent calculation by Hood et aL{12), other key quantities in DFT, namely 
the exchange-correlation hole Uxc and the exchange-correlation energy density 
Cxc, have not been investigated with Monte Carlo methods. These quantities, 
however, are important in understanding the success of the LDA beyond its 
formal limits of validity, and play a key role in constructing more accurate 
approximations to E^c- A better knowledge of these quantities is therefore 
crucial for a better understanding of the performance of the LDA and various 
corrections to it such as GGAs, and can guide the construction of more accu- 
rate functionals. Unlike V^c which can be directly obtained from the electron 
density (by inversion of the Kohn-Sham equations ([T0| , |^)) evaluating e^c and 
Uxc is more demanding. These quantities are derived from an adiabatic con- 
nection procedure in which one scales the Coulomb interaction by a factor A 



while keeping the density fixed at the ground state density of the system un- 
der consideration. To extract Uxc and exc from Monte Carlo data one therfore 
needs to calculate not only the ground state many-body wavefunction of the 
fully interacting system (A = 1), but also the many-body wavefunction in the 
range < A < 1. 

Within variational quantum Monte Carlo, we have developed a new scheme 
for realizing the above adiabatic connection procedure which allows us to ex- 
tract rixc and e^c from Monte Carlo data. Our method is based on a constrained 
optimization of the many-body wavefunction at different Coulomb coupling 
constants using the technique of variance minimization p^) . In this paper 
we will discuss aspects of our method and illustrate it with a first application 
to an electron gas exposed to a cosine-wave potential. For this system we 
calculate the exchange-correlation energy, exchange-correlation energy density 
and exchange-correlation hole, and compare our findings with those obtained 
from the LDA and the most commonly used version of GGA (^, |15|). 



2 The Adiabatic Connection 

The idea of an adiabatic connection to determine E^c has been developed by 
several authors (16-18). Here we closely follow the review by Parr and Yang 
d^). We consider a system of interacting electrons in the presence of an 
external potential Kx(r) and characterized by the Hamiltonian (atomic units 
are used throughout, with e = Ti = m = 1) 

H = f + V,e + Vex (1) 

with 

^ 1 

T = T.~'^' (2) 
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Ve. = Y.VeM) (4) 
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In the Kohn-Sham formulation of DFT the problem of finding the ground state 
energy of this system is exactly mapped onto one of finding the electron density 
which minimizes the total energy functional 



E[n{v)] = To[n(r)] + EhW)] + j dvVex{v)n{v) + E,,[n(r)] (5) 



Here Tq is the kinetic energy of a fictitious non-interacting system of elec- 
trons having the same electron density n(r) as the interacting system and 
EhIu] is the Hartree (electrostatic) energy. The exchange-correlation energy 
functional -E^ci^] is usually defined by equation (5) and contains all the many- 
body terms not considered elsewhere in (5). 

An exact expression for E^c is obtained by scaling the electron-electron 
interaction with a factor A and varying A between 1 (real system) and (non- 
interacting system). The exchange-correlation functional Ercc is then given 
by (i 

E.c[n] = C d\< ^^iKel*^ > -EhH (6) 



where is the anti-symmetric many-body wavefunction which minimizes 

F^=<t + AV;e> (7) 
under the fixed-density constraint 

< ^^|n(r)|^^ >= n(r) (8) 
and 72 (r) is the density operator 

N 

n{r)=Y.5{v-Vi). (9) 

i=l 

A minimum for F^ always exists (|20|) and, except under some unusual condi- 
tions (|2lD, \1'^ can be obtained from the following Schrodinger equation 

[f + XVee + V^]^^ = Hx^x = E^^^ (10) 

with 

N 

v' = Y.y\^i)- (11) 

i=l 

The potential V^^(r) at point r is a Lagrange multiplier corresponding to the 
fixed-density constraint at that point. As A varies between and 1, l^^(r) must 
be adjusted such that the electron density remains fixed at nir). At A = 1, 
coincides with the actual external potential Kx(r) while at A = 0, it coincides 
with the Kohn-Sham effective potential, 

V''=\v) = V,ff{v) = Ve.{v) + VH{r) + F,,(r) (12) 

where Vh is the Hartree (electrostatic) potential 



(13) 



and 

is the Kohn-Sham exchange-correlation potentiaL Note also that corre- 
sponds to the Slater determinant of the exact Kohn-Sham orbitals correspond- 
ing to the density niv). 

The adiabatic expression (6) allows us to obtain several useful decomposi- 
tions of Exc- Inserting (3) in (6) gives 

E^c\n{v)\ = - J dr J dr (15) 

where n^c is the density-functional exchange-correlation hole defined by (|^) 

n(r, r') = n{r)n{r') + n{r)nxc{r, r') (16) 

Here n(r, r') is the diagonal part of the two-particle density matrix averaged 
over A, 



') = [^dXn\r,r') (17) 
Jo 



n(r, r 
and 
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n\r, O =< v&^l E E ^(r - - r,)\^' > . (18) 

Integrating (15) over r' yields 

E,,[n{r)] = I rfre,,(n[r],r) (19) 

where e^c is the exchange-correlation energy density derived from the adiabatic 
connection procedure 

e.,([n(r)],r)= /' rfA ei([n(r)], r) (20) 
Jo 

with e^(, given by 

x /r / M X T 1 1 (^(r - Tj) , ^ 1 /■ , ,n(r)n(r') , ^ , 

eLiHr)],v) =< ^aI-EE^— TtI^a > '2 7^' r - r1 ^^^^ 

i=l jj^i I i I II 

For further reference we note that corresponds to the density func- 

tional exchange hole n^. The corresponding excahange energy density is = 
e;^^° and the correlation energy density is given by Cc = e^c — e^. Note, how- 
ever, that the excahnge-correlation energy density, and hence its exchange and 



correlation components, are not uniquely defined quantities since we can al- 
ways add to txc any function which integartes to zero without affecting the 
exchange-correlation energy. Our definition of these quantities emerges in a 
natural way from the adiabatic connection. An alternative definition of the 
correlation energy density has been suggested by Baerends and Gritsenko (p^) 
and by Huang and Umrigar (p3[). 



3 Quantum Monte Carlo realization 

Given an interacting many-body system with ground-state density n(r), the 
main ingredient for evaluating n^c and txc is the many-body wavefunction 
\E'^ for a number of systems corresponding to different values of the coupling 
constant A satisfying the fixed-density constraint. In this section we describe 
our variational quantum Monte Carlo algorithm for obtaining "^x. 

3.1 Variational Monte Carlo 

In variational Monte Carlo calculations (0) one starts off with an explicit 
parameterized Ansatz for the ground-state many-body wavefunction of the 
system under consideration. The total energy of the system is then calcu- 
lated as the expectation value of the Hamiltonian H with respect to the vari- 
ational wavefunction Monte Carlo integration is used to perform the 
multi-dimensional integrals required for evaluating this expectation value and 
the variational parameters in are adjusted until an optimized wavefunc- 
tion is obtained. The state-of-the-art method for performing the optimization 
procedure is the variance minimization scheme (|T^, |T^. In this scheme one 
minimizes the variance of the local energy H'^t/^t (rather than expectation 
value of H) with respect to variational parameters over a set of particle config- 
urations. The use of energy optimized wavefunctions may give unsatisfactory 
results when quantities other than the energy are evaluated, while minimiza- 
tion of the variance tends to give a better fit for the wavefunction as a whole, 
so that satisfactory results are obtained for a range of quantities including both 
energy and electron density. The electron density plays a central role in the 
adiabatic connection procedure making variance minimization a more suitable 
choice for optimizing 

3.2 Fixed-density variance minimization 

We consider an N-electron system having ground-state density n(r). At a given 
coupling constant A the corresponding many-body wavefunction satisfies 



equation (10). Therefore, at an arbitrary point R = (ri, r2, . . . r^) in the 3A^ 
dimensional configuration space of electron coordinates, we have 



i7^^''(R) 



- E^ = 0. (22) 



In conventional variance minimization calculations (|T^) (i.e. the unconstrained 
A = 1 case), the above property is used to find an overall fit to "if (we drop the A 
superscript for simplicity). The procedure is to determine the parameters {a} 
in the trial function \E't(R, {ct}) by minimizing the variance of local energy 



a^= I dR 



I^t(R)P (23) 



where E[iI)t] is the expectation value of the Hamiltonian. 

The above unconstrained optimization cannot be directly applied at inter- 
mediate values of A for which the Hamiltonian contains the unknown potential 
Vx- We found, however, that a simultaneous determination of \E'a and Vx can 
be achieved by performing the following constrained optimization. We assume 
that the trial many-body wavefunction results in the electron density n^(r) 
and expand both n'^(r) and the ground state density n(r) in a complete and 
orthonormal set of basis functions {fs} 

n{r) = J2nJs{r) (24) 

s=l 

n\r) = j:<fs{r) (25) 

s=l 

where A''^ is a cut-off chosen such that the above expansions converge to n{r) 
and ^^(r) within a specified accuracy. Subsequently, we define the modified 
penalty function /x^ 

^^ = a' + Wj2[n,-n^] (26) 

s=l 

where is a weight factor the magnitude of which determines the emphasis 
laid on the fixed-density constraint. The above penalty function reaches its 
lower bound (of zero) if and only if \E'^ is the exact many-body wavefunction 
satisfying the fixed density constraint (within the accuracy set by Nd) and is 
the corresponding exact potential. Hence minimization of /x^ will, in principle, 
result in the simultaneous determination of \E'^ and V^. In practice, however, 
our constrained search is restricted to a sub-space of many-body wavefunctions 
and minimization of /i^ yields an optimal fit to \E'^ and a corresponding optimal 



fit to V^, the deviations of wliich from the exact reflect the errors in the 
many-body wavefunction. 

Our numerical implementation of the above scheme works as follows. We 
start off with an initial guesses \E'o for the many-body wavefunction and a 
corresponding guess for V''^. A fixed number Nc of statistically independent 
configurations Rj are then sampled from |\E'oP and the Monte Carlo estimator 
of over these configurations is evaluated 
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^m^ = J2{El{R,)-<El>) 



UJi 
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< El > the average energy 

<El >=J2EL{Ki 
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The expansion coefficients of the electron density, n^, are evaluated from 



Nc N 
i=l k=l 



UJi 



Nc 



(27) 

(28) 
(29) 

(30) 
(31) 



where denotes the coordinates of the electron k belonging to configuration i. 
Finally, we vary the parameters in \E'^ and V''^, using a standard NAG routine 
for optimization, until /i^ is minimized. We found that setting W equal to 
the number of configurations results in a satisfactory minimization of both the 
variance in energy and the error in electron density. 

Following (|T^) we set the reweighting factors uJi in equations (27) and (30) 
equal to unity in order to avoid a numerical instability in the variance mini- 
mization procedure which occurs for systems with a large number of electrons 
(these factors, however, are included in calculating the expansion coefficients 
of the electron density). The above fixed-density variance minimization is then 
repeated several times until the procedure converges. 



4 Cosine- wave jellium 

We performed adiabatic connection calculations for the inhomogeneous spin- 
unpolarized electron gas with average electron density rio = 3/(47rr^) corre- 
sponding to = 2. In the QMC simulations we model this system by a finite 



system of = 64 electrons satisfying periodic boundary conditions in a FCC 
simulation cell. Density modulations can be induced by applying an external 
potential of the form cos(q.r) where, because of periodic boundary condi- 
tions, q is restricted to be a reciprocal lattice vector of the simulation cell. 
Alternatively, we can fix the ground-state electron density a priori and use 
our fixed-density variance minimization method to obtain the corresponding 
many-body wavefunction at a given Coulomb coupling constant which pro- 
duces the specified density. In the calculations reported here we chose this 
second option, with the "target" electron density for the system generated in 
the following way. We expose the non-interacting electrons (i.e. the A = 
system) to the potential ^(r) 

y(r) =T/qCOs(q.r) (32) 

with Vq = 2.084:6% and q = 2B3 . Here e% is the Fermi energy of the un- 
perturbed electron gas, B3 is a primitive vector of the reciprocal (simulation) 
cell with I2B3I = l.llkp, and kp is the Fermi wavevector. We then solve the 
following self-consistent single-particle Schrodinger equations 

[~V^ + Veff](l>i = eicj>i (33) 

with 

K//(r) = V{v) + Vnir) + V^\^\v) (34) 
to obtain the electron density 

N/2 

n{v)=2Y.\m? (35) 

We define this density to be the exact ground-state density of our interacting 
system. In this way, the single-particle orbitals 0j are by construction the ex- 
act Kohn-Sham orbitals and their Slater determinant corresponds exactly to 
the many-body wavefunction at A = 0. Having obtained this non-interacting 
i;— represent able density we then perform fixed-density variance minimization 
to produce variational many-body wavefunctions at non-zero Coulomb cou- 
pling constants (including the ground-state many-body wavefunction at A = 1) 
which reproduce this density and (variationally) satisfy the Schrodinger equa- 
tion (10). Once the \E'as are obtained, we use the Monte Carlo Metropolis 
algorithm to evaluate the required expectation values and perform a numeri- 
cal coupling constant integration using Gaussian quadrature. 



4.1 Many-body wavefunction 



The quality of a variational quantum Monte Carlo calculation is determined by 
the choice of the many-body wavefunction. The many-body wavefunction we 
use is of the parameterized Slater- Jastrow type which has been shown to yield 
accurate results both for the homogeneous electron gas and for solid silicon 
([T^ ) (In the case of silicon, for example, 85% of the fixed-node diffusion Monte 
Carlo correlation energy is recovered). At a given coupling A, is written as 



D^D^ exp 



1>J 



(36) 



where = Ir — r^ l and and are Slater determinants of spin-up and spin- 
down Kohn-Sham orbitals respectively, is the two-body term correlating 
the motion of pairs of electrons and cxj denotes the spin of electron i. Finally, 
X'*' is a one-body function which is absent in the homogeneous electron gas but 
is crucial for a satisfactory description of systems with inhomogeneity. Both 
and contain variational parameters. We write m'^ as (|I4D 



u\r)=4{r) + f\r), 
where Uq is a fixed function and /'^ is given by 



(37) 



fir) 



B 



2 +r){Lws-ry + r\Lws-ryEfio(ytTi(f) 0<r<Lws 
r > Lws 



(38) 

where and are variational coefficients, Ti is the Ith Chebyshev polyno- 
mial, and 



2r-L 



ws 



(39) 



In the last two equations Lws is the radius of the sphere touching the Wigner- 
Seitz cell of the simulation cell. 

The fixed part of at full coupling constant A = 1 is the short-ranged 
Yukawa form (|T^) 



(40) 



where is fixed by the plasma frequency of the unperturbed electron gas 



A' 



(41) 



and is fixed by imposing the cusp condition (|^) leading to F^_^^_. = y (2^4^) 

for parallel spins and F^. = \J\a^) for anti-parallel spins. Lq is a cut-off pa- 
rameter chosen so that Uq{Lws) is effectively zero and is set equal to 0.25Lvys 
in the present calculations. In the case of the unperturbed electron gas, scaling 
arguments (^) applied to the Hamiltonian (10) result in the following relation 
for the exact many-body wavefunction at coupling constant A 

^,"^(ri, r2, . . . r„) = C"^^r^(Ari, Ar^, . . . Ar„) (42) 

where '^'tr^ is the ground-state wavefunction of a homogeneous electron gas 

' s 

with the density parameter = Ar^ and is a normalization constant. For 
the unperturbed electron gas (x = 0) imposing condition (42) on the fixed-part 
of our Slater- Jastrow wavefunction yields 




where = X^^'^A^, F-^ = A~^/^F^. We note that with the above choice for 
A'^ and F^ the A— dependent cusp conditions are automatically satisfied. The 
electron density is modulated only in the B3 direction and hence both the 
one-body part of the Jastrow factor and can be expanded as 

M 

x\r) = J2 X^rnBs) cos(mB3.r) (44) 

m=l 
M 

V\r) = ^ V^mBs) cos(mB3.r) (45) 

m=l 

The electron density is expanded in a similar way (with the inclusion of the 
m = term) . We use 7 Fourier coefficients in the expansion of electron density, 
6 Fourier coefficients in the expansions of xx and Vx (only the first four coeffi- 
cients turned out to be significantly different from zero), and 8 coefficients (for 
each of the spin-parallel and spin-antiparallel cases) in the two-body term. 

5 Results and discussion 

We performed adiabatic connection calculations for cosine-wave jellium us- 
ing six values of A: 0,0.2,0.4,0.6,0.8,1. The many-body wavefunctions for 
A > were optimized by fixed-density variance minimization using 10000 in- 
dependent A^— electron configurations at each A. These configurations were 
regenerated several times. The weight factor in expression (27) was set equal 




Figure 1: (a) Electron density for different values of A plotted along ( (the direction 
of inhomogeneity) . (b) The A— dependent spherically averaged exchange-correlation 
hole for an electron sitting at ^ = 10.85 a.u. 

to 10000 in order to obtain a satisfactory minimization of both the variance in 
energy and the error in electron density. Once was optimized, quantities of 
interest were accumulated with the Metropolis Monte Carlo algorithm using 
500000 statistically uncorrelated configurations. 

We found that our method results in electron densities n^(r) which deviate 
from the reference density by less than 1%. This is shown in figure 1(a) where 
the density is plotted as a function of A along a line parallel to the direction in 
which the external potential varies (we call this the ( direction). While the den- 
sity is fixed, all other physical quantities vary smoothly and monotically with 
A. As an example we consider the spherically-averaged exchange-correlation 
hole 

n^v, s) = ^ f dr'n^r, r'), Q : |r - r'| = s (46) 

An Jn 

as a function of A. In figure 1(b) this quantity is shown around an electron 
sitting at one of the maxima of the electron density {( — 10.85 a.u.). The 
A = curve corresponds to the spherically-averaged exchange hole. The ex- 



change hole is relatively shallow and negative everywhere. As the interaction 
is switched on, the hole around the electron becomes gradually deeper. The 
spherically-averaged hole obeys the sum-rule 



47r / s'^h^{r,s) = -1 (47) 



and the deepening of the hole for A > is compensated by the fact that the 
hole becomes slightly positive far away from the electron. Note that the hole 
does not "narrow" as it deepens, but actually broadens. In evaluating p^^ we 
expanded the exchange-correlation hole in a double Fourier series, sampled the 
corresponding expansion coefficients and subsequently performed the spherical 
averaging. Our calculated p'^^ does not satisfy the Kimball cusp-condition 



291) because of the finite number of plane-waves in its Fourier expansion. As 
a result p'^^ has zero slope at s = 0. We note, however, that this deficiency 
does not affect -E^^c and e;^^ because these quantities are evaluated directly from 
equations (6) and (21). 

Before discussing our findings for this system, we would like to pause and 
give a short outline of the errors present in our simulations. First of all, the 
small (< 1%) deviations of the electron density at different A from the refer- 
ence density ?7.(r) will induce errors in the adiabatically calculated quantities 
such as the correlation energy density. By recalculating the exchange energy 
density with a density which deviates from the reference density by 1% and 
extrapolating the resulting deviation e2^[n(r)] — ex[n{r) +6n{r)] to the correla- 
tion energy density, we estimate the errors in Cc due to these density deviations 
to be also ~ 1%. 

Further, there are two other kind of errors in our calculations: (i) statistical 
errors; and (ii) finite size errors (i.e. those caused by the fact that we are 
using a finite number of electrons to model a supposedly infinite system). 
With 500000 configurations used in sampling all physical quantities, we found 
statistical errors to be unimportant, except for the exchange-correlation hole. 
By evaluating the exchange hole both directly and by Monte Carlo sampling, 
and assuming that the errors in rixc for A 7^ are similar, we estimate the 
statistical error in this quantity to be less than 7%. Another source of errors 
is finite size effects. These errors occur because a finite simulation cell is used 
to model an infinite system, with the Coulomb interaction energy evaluated 
using the Ewald formula (|26|). The use of a finite simulation cell with periodic 
boundary conditions affects the wavefunction, of course, and the use of Ewald 
interaction also produces a Coulomb finite size error in the interaction energy 



([2B|, P?]). We found the effect of the finite cell on the exchange-correlation hole 



to be unimportant, except for the asymptotic behavior of this quantity which 
cannot be correctly described with the present system size. Coulomb finite size 
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Figure 2: (a) VMC (solid lines) and LDA (dashed lines) snxc{T^,s) plotted for an 
electron moving along C. Arrows on the electron density (plotted on top), mark the 
position of the electron, (b) Exact shx{Tc,s) plotted in the same direction and at 
the same points as in (a) (solid lines) and the corresponding LDA approximation 
(dashed lines). 



effects do not significantly affect the many-body wavefunctions, and hence their 
effect on quantities such as the electron density and the exchange-correlation 
hole is negligible. The exchange-correlation energy density, however, is directly 
affected by Coulomb finite size effects since in evaluating e^^ the l/|r — r'| 
Coulomb interaction in (21) is replaced by the periodic Ewald interaction. By 
calculating the exchange energy density of the homogeneous electron gas using 
our finite simulation cell and comparing it with the exact result (|28D 

0.45805 , 
-no (48) 



we estimate the total finite-size error in to be of the order of —2 x 10 ^ 
a.u; the errors in the correlation energy density Cc is expected to be somewhat 
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Fi gure 3: Exchange (solid line) and correlation (dashed line) contributions to e.xc 
plotted along C^. 

smaller. 

We now turn to our results for n^f. and e^c- The spherically averaged 
exchange-correlation hole, fi;^^(r, s), obtained from our adiabatic calculations 
is shown in figure 2(a) together with the LDA approximation ( pSj) to this 
quantity 

n^f ^(r, s) = n(r)(/°-(n(r), s) - 1) (49) 

where g^""^ is the A-averaged pair-correlation function of a homogeneous elec- 
tron gas with density n(r) (we use the Perdew-Wang parameterization of g 
(|25|)). In this figure we plot nxdi^, s) for an electron moving along (, and hence 
fully experiencing the strong variations in electron density. The hole is shown 
multiplied by s so that the area under each curve is directly proportional to the 
exchange-correlation energy per electron e^c/n. At ( = 1.30 a.u. the electron 
density is very low and rixc is shallow. As electron moves to higher densities 
(C = 5.6 and ( = 10.85 a.u.) the hole becomes deeper and its asymptotic tail 
less pronounced. Unlike the LDA hole which depend only on the local density 
n(r) and is dug out of a homogeneous electron gas of that density, the VMC 



hole depends on the density everywhere in the vicinity of r. At C = 1-30 a.u, 
where the electron density is very low, the LDA "probes" only this density and 
for this reason the LDA exchange-correlation hole is very different from the 
VMC hole, even close to the electron. As the electron moves to the high den- 
sity region, the LDA description becomes more satisfactory and at ( = 10.85 
a.u., where the electron density has a maximum, the agreement between the 
LDA and the VMC hole is rather good. In figure 2(b) we compare the ex- 
act exchange hole (obtained from our exact Kohn-Sham orbitals) for the same 
electron positions with the LDA hole given by 



Jl(fcir(r)g) 



(50) 



where fc_F(r) = (37r^n(r))^/^ is the local Fermi wavevector and ji the first order 
spherical Bessel function. Once again, the LDA description is unsatisfactory 
at low densities but improves as we move to the high density region. 

Next we consider exchange-correlation energy densities. In figure 3 the 
exchange and correlation contributions to this quantity are shown. The differ- 
ences ejf^*" — e^'^^ and e^^*" — e^^^ are shown in figure 4 (a). The difference 
in Cc follows the variations in electron density and is largest at points where 
n(r) has a maximum. The differences in Cx shows a more complicated structure 
and ejf*^'" < e^^"^ everywhere in the system. This result is in line with the 
well-known fact that LDA almost always underestimates the exchange energy 
of an inhomogeneous system. However, because of the finite size errors, we 
expect the true to be slightly (~ 2 x 10~4 a.u.) less negative than ejf^'" so 
that < G^^^ must holds in most points of the structure but not necessarily 
everywhere. 

The GGAs for exchange and correlation of a spin-unpolarized system are 
written as ( |T5|) 

Ex = jdv n{v)eT'f{n{v))Fx{s) (51) 

E^ = J dr n{r) [e™*-^(n(r)) + H{n{r),t)] (52) 

In the above equations e""*-^ and e™*/ are the exchange and correlation energies 
per particle of a uniform electron gas with density ri(r), t = | Vn|/2/cs(r), s = 
\Vn\/2kp{r) with ks the local Thomas- Fermi wavevector ks{r) = ^4:kp{r)/n. 
In analogy with (19) one may define the GGA and Cc as 

e^^^{r)=n{r)erf{n{r))Fx{s) (53) 

and 

ef«^(r) = n(r)[er/(n(r)) + H{n{r),t)] (54) 
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We note that the above quantities do not directly correspond to the physi- 
cal and Cc as defined by equation (21) (i.e. via the exchange-correlation 
hole). Nevertheless, they present pointwise corrections to the LDA exchange 
and correlation energy densities and for this reason we found it interesting to 
compare them with our VMC results. In figure 4(b) we show e^*^*" — e^*-^"^ 
and e^^*" — e^*^"^. The GGA results were obtained from the ground state den- 
sity n(r) using the Perdew-Burke-Ernzerhof scheme (|T5|) , which we found to 
give results slightly different from PW91 (^. Note that the difference between 
the VMC and the GGA exchange energy densities is significantly smaller than 
that between the VMC and the LDA exchange energy densities, indicating 
that the GGA improves upon the LDA in describing this quantity. As for Cc, 
the differences are of the same order as for the LDA, although the shape is 



different. It is interesting that the LDA errors in Cx and Cc partially cancel 
each other, even on a local scale, but that these cancellations do not occur for 
the GGA. In summary, the GGA seems to do a good job in improving the LDA 
description of the exchange energy density but is less successful in the case of 
the correlation energy density. The resulting exchange correlation energy per 
electron obtained from VMC, LDA and GGA are E}^f^/N = -0.328 ± 0.009, 
E^c^/N = -0.3296, E'^^^/N = -0.3347 a.u. 
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